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Abstract—The problem of diffraction of an external electro- 
magnetic field by a locally inhomogeneous body placed in a rect- 
angular waveguide with perfectly conducting walls is considered. 
The formulated problem is reduced to a volume singular integral 
equation. The problem is solved with the use of the numerical 
collocation method. Volume singular integral equation method 
is applied to analyze the problem for structures of a complex 
geometric shape. Since the problem is highly computer- 
intensive, it is solved with the use of parallel algorithms realized 
on a supercomputer. 


Index Terms - Maxwell equations, volume singular integral 
equation, Galerkin method, numerical results. 


I. INTRODUCTION 


The study of problems of diffraction of electromagnetic 
waves in waveguides leads to the solution of three-dimensional 
vector problems to the complete electrodynamic formulation. 

At present, solution of such problems is one of the top- 
ical tasks in electromagnetics. The solution that is based on 
mathematical methods and guarantees the accuracy suitable for 
practice at the electrodynamic level of rigor necessitates a very 
large amount of computation, which cannot be implemented 
even on the state of the art supercomputers. Solution of inverse 
electrodynamic problems for a complex system of surfaces and 
bodies in the resonance frequency region is the most urgent 
task. This problem arises when the parameters of nanocom- 
posite materials and nanostructures are to be determined. 

Solving the aforementioned problems with the help of 
expensive applied software packages using traditional finite 
difference techniques or finite element methods does not 
always lead to sufficiently accurate results. 

Since exact solutions of diffraction problems can be ob- 
tained for a limited set of bodies with a regular geometry, 
it is important for many applications to develop various 
approximate numerical methods valid for bodies of arbitrary 
shapes. Thus, it is necessary to develop new methods for 
solving this kind of problems. 

One of the promising methods is the method of volume 
singular integral equations (VSIE), by means of which the 


boundary problem is reduced to solving a VSIE [1], [2], [3], 
[4]. Generally, the obtained integral equation can be solved 
only by numerical methods but, since the dimensionality of 
the problem is reduced due to the reduction of the problem to 
a surface integral, numerical calculation is significantly sim- 
plified. Solving such problems with an accuracy appropriate 
for practice requires a considerable amount of computational 
effort. 

At present, a rapid progress in the computer technology 
facilitates wide application of computer simulation methods 
for solving such problems for screens of canonical shapes. 

In this work, the development of a software based on 
subhierarchical method is proposed and implemented. The 
method presented in this paper makes it possible to solve the 
problems for inhomogeneous bodies of complex geometric 
shapes using the results of the solution of the problem for 
bodies of basic shapes. 


II. FORMULATION OF THE PROBLEM 


Consider the following diffraction problem. In Cartesian 
coordinates, let P = {x : 0 < a4 < a0 < a < 
b, =œ < z3 < + co} be a waveguide with perfectly 
conducting surface OP (Fig. 1). The waveguide contains 
volume body Q (Q C P being a region) characterized by 
constant permeability uo and positive 3 x 3 tensor permittivity 
function (x). The components of ê(x) are bounded functions 
in region Q, inverse tensor ê € L,,(Q) exists in Q, and its 
components are also bounded in Q . 


Fig 1. The body @ in the waiveguide P. 


The region ðQ is assumed to have a piecewise smooth 
boundary Q . We assume, in addition, that the body Q does not 
touch the waveguide walls, i.e. OQ N OP = Ø. Beyond P\Q, 
the medium is isotropic, homogeneous, and characterized by 
constants €9(> 0), uo(> 0). 

It is necessary to determine the electromagnetic field’s 
E,H € Lə1oc(P) (and, hence, FE, H € Lz (Q)) excited in 
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the waveguide by an external field with, the time dependence 
exp(—iwt), where w is the circular frequency. 

In region P C R, the standard differential operators 
grad, div and curl are understood in the sense of generalized 
functions. 

We seek weak (generalized) solutions to the system of 
Maxwell equations 


curlH = —iwéEH 


curl = iwuoH, x € P. ` (1) 


These solutions must satisfy the conditions at infinity: fields 
E and H for |x3| > C and rather large C > 0 are represented 


as 
E 
(a)=( æ )+ 
“ i ie =v 
+ Reexplinpleal) ( p Ipes — i% V2 2 
Sepa Vay) x €3 


7 —iweg(Ve2Il,) x C3 
+ erin? |x ( 2 
2° pliyplzal) S — i Vorhy (2) 


(“+” and “— correspond to +co and —oo, respectively), 
where ĉi 2 3are units vectors in the Cartesian coordinate frame, 


IP = yk? — AY, Img? > 0 or Imo? = 0, kor? > 0 


and [ my! i) = 0; ko’ ) > 0 sets of eigenvalues and eigen- 
functions (orthonormal in Lə (II)) of 2D Laplace operator — A 
in the rectangle I] := {(21,2%2):0< 2, <a,0< z2 < b} 
with the Dirichlet and Neumann conditions, respectively; 
ke = €9 Uw. 

The coefficients RS and QEP appearing in the right hand 
side of (2) are estimated to have the behavior 


E® 
H? 


RS), QGP = O (p™) ,p > œ (3) 


for certain m € N. 

E? and H°, appearing in (2) are known incident fields, and 
they are the solutions of the boundary value problem in the 
absence of an inhomogeneous body Q. 

From a physical viewpoint, conditions expressed in (2) 
that the scattered field is the superposition of normal waves 
emanating from the body. Conditions (3) provide for the 
exponential convergence of series (2) and the possibility of 


their term-wise differentiation with respect to x; performed 
an arbitrary number of times. 

The field & must satisfy the boundary condition on the 
waveguide’s walls 


(4) 


Relationships (1)-(4) for the field & yield the integrodiffer- 
ential equation [4] 


E,|ap = 0. 


E (x) = E’ (x) +b | Ge (r) (x i E (y) dy+ 
Q 
graddiv | Gp (r) (Ha — i E (y)dy EQ (5) 


Q 


where T is the unit tensor. i 
The components of the diagonal Green’s tensor Gz = 
diag (Gz, G}, G3) have the following form [4], [6], [7], [8]: 


exp ( —Ynm [£3 — y3|) 
a ao 
z uno 


m 
b Y2, 


Tn nM Tn 
cos — %1 sin ——2%2 cos — yı sin 
a b a 


D3 > exp ( —Ynm [£3 — y3|) . 
n=1m=0 Yam (1 + dom ) 
nm 


` Tmn Tnm Tmn 
sin — £1 COS — T2 sin — Yı COS —— Y2, 
a b a b 


25 5 exp ( —Ynm a ysl) 


Becerra 


awn i. 7m . mm . TM 
SIN — T1 SIN —— T2 SN — Yı SIN —— Y2. 
a b a b 
i 2 2 
In these expressions, Yam = \/ (=) + (==) — kê the 
branch of the square root is chosen such that IMYynm => 0. 
Equation (5) can be solved by means of various numerical 
methods. In this study, we apply the collocation method, 
because the use of the Galerkin method leads to more cum- 
bersome formulas and calculations. 


IHI. THE COLLOCATION METHOD 


For the equation Ad = f,(¢, f € X), with linear bounded 
operator A : X — X in Hilbert space X, we consider the 
collocation method that is formulated as follows. Approximate 
solution ¢, € Xn is determined from the equation P, Ad, = 
P,f. Here, n E Xn (Xn is an n-dimensional subspace of 
space X), and Pa : X — Xn 1s result of the projection of 
X onto a finite dimensional subspace, by using the operator 
defined below. 

Let us split the region Q into elementary subdomains 
Q, with piecewise smooth boundaries 0@; so that the condi- 
tions Q; N Q; = 0 are fulfilled when i 4 j andQ = UQ;. 


7 
In each subdomain Q; , we choose the collocation point 
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l«reQd; 
0,2 Qs 


Let subspaces be linear hulls of basis functions: Xn 
span{v, ..., Un}. 

We require that the chosen basis functions satisfy the 
approximation condition 


(node) x’. Consider the basis functionsv; = 


vVreX lim inf |e — z|| = 0. 


Let us define operator Pp X — Xn as follows: 


(Pho) (z) = olz’), £ E Qi 

Note that, in this situation, the values of functions (Pao) (x) 
for x € OQ; are not defined, but this is not important, because 
X = Lə in our Case. 

The equation P, Ad, = Pa f is equivalent to the following: 


Abie Sfe% 7H hm: 
Let us represent the approximate solution as a linear com- 
n 
bination of basis functions: 6, = >) CkUk 


k=1 
The substitution of this representation into the scheme of 
the collocation method yields the system of linear algebraic 
equations for unknown coefficients cg: 


ie 4 ere 


S cx(Avg) (a!) = f(a! 
k=1 


Now, let us construct a scheme for solving the integral 
equation by means of the collocation method. 
We formulate the method for integrodifferential equation 


(5). Assume that the tensor (42 sd ) 
= 
E Doo (Q). 


Let us introduce the notation 


A —1 
é- (2-7) ws 
EQ 


and pass to the next equation 


is invertible in Q, 


SORES, 
EO 
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AJ =Ê) Jæ) -rè | Gee») J(u)dy- 


graddiv | Cred =J (z) EQ. (6) 
Q 


Equation (6) can be represented as a system of three scalar 
equations: 


S &(e) = 6B | ele. v) y)dy- 


Q 


O ‘ 
sotive | Ĝele y) IU) = E” (g), L=1,2,3. (7) 
l 
Q 
We define the components of the approximate solution Jp = 
(Ji, JŽ, J?) as follows: 


J = Oki. aed, = 0. 
k=1 k=1 


k=1 


where f} are the basis step functions. 
Below, we construct functions f$. 
parallelepiped: 


Assume that @ is a 


Q = {x : a1 < T1 < Q2, bi < £2 < bg, Cy Ea < Cz}. 


Let us split Q into elementary parallelepipeds: 


Wie oie = 2h Spit a 


L3.m < T3 < DS obi f 


—a bə — b 
2 Ly piei 2 1 








Ti k = Q1 + l, 





where k,l, m = 0,... n — 1. 
We obtain the following formulas forf}),,, i = 1,2,3: 
1, x € kim 


a = = 
Tri i 0, r g Ikim 

The constructed set of basis functions satisfies the necessary 
approximation condition in L? = Lə x Le x In. 

In order to construct matrix elements, we integrate the 
components of the Green’s tensor over the parallelepiped 
es os { (21, £2, £3) : t1 =< ~ < t1 eed 19 < r < 
i2 + lig < r < i3 + 1} and denote these components 





G+, G2 G?. Then, we obtain 
C > > J Jam (23) (z3) cos (n.X 1) sin (mXə2) sin Rii 
L 72 Pee k y2 nnm f i 2 





mH 
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2H H 
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H 
x sin (nH; (i; + 0.5)) cos (mH (i2 + 0.5)) sin (= :) P 





2H H 
2Hı Po ae li (nXı)sin (nH; (i, + 0.5)) sin (=) 


H 
sin (nXı)sin (mXə)sin (=) 


J fige 


rome em i 





H 
x sin (nH; (i1 + 0.5)) sin (mHə (i2 + 0.5)) sin (= z) 


Here, 
ee, eg eee ye 
a b a 
h h 
pe lias mao 
b a b 
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tı = ili, z2 = johe, yi = tihi, Y2 = təhħ2, 


wT 


a 


exp(—(x3 — (i3 + 1)h3)Ynm) 
—exp(— (z3 a igh3)Ynm), 
E (ee his 


exp(—(izhs a £3) Ynm) 
—exp(—((is + 1)h3 — £3)Ynm); 
gzz < ishs 


LAs) = 


2 — exp(—(x3 — i3hs)Ynm) 
—exp(—((%3 + 1)h3 — £3)Ynm); 
i3h3 < £3 < (i3 es 1)hg 


The values of the integrated components of the Green’s 
tensor at a collocation point are obtained by summating 
slowly convergent series and their second order derivatives. 
The computation of slowly convergent series is accelerated by 
separating the singularity. 

It is convenient to represent the extended matrix obtained 
according to the collocation method in a block form: 


Ay. Aig Ag | Bi 
A2 Áz A233 | Bo 
A31 A32 A33 | B3 


Elements Axı and Bk are determined from the relationships 


R (az) — Syah? J G*(a,,y) fily)ay 
Q 


a { 
-a a a 5 Glez y) fi U)dy 8) 


By, = Ey (2:) (9) 


where the coordinates of the collocation points are 
Ti = (Zil, 072,053) , Zil = (i1 + 0.5) hy, 
2 = (i2 + 0.5) ho, ziz = (13 + 0.5) h3, 

ble 128s U1 as a ig 3 Ua = L 
Thus, we have obtained the necessary formulas for calculat- 
ing the matrix coefficients of the collocation method applied 


to solve the volume singular integral equation in the problem 
of determination of the permittivity of a material. 


IV. THE SUBHIERARCHIAL ALGORITHM 


Let Q = 4a + ay < zti < @, bi < a < 
bo, C1 < z3 < co} be a filled waveguide section 
(a, = 0, ag = a, bı = 0, bə = b, Cı = 0, = Cc). The al- 
gorithm of calculation of the electromagnetic field inside a 
filled waveguide section has been described in the foregoing. 
Let us consider the algorithm of calculation of the electromag- 
netic field in a rectangular waveguide for a body of a complex 
geometric shape. We use the matrix obtained according to 
the collocation method for a filled section. The problem of 
diffraction by a body of a complex shape can be solved when 
the body entirely fits into the considered waveguide section 
and consists of grid elements [5], [9], [10], [11], [12]. 

The subhierarchical method makes it possible to compose 
a matrix for determination of the electromagnetic fields inside 
a body of a complex shape with the use of the matrix calcu- 
lated for the complete waveguide section. In the constructed 
structure, we introduce a new numeration of elementary paral- 
lelepipeds belonging to the structure of a complex shape and 
form a new grid. This grid will be employed for calculating 
the field on the body of a complex shape. 

We solve the system of linear algebraic equations for the 
matrix composed with the use of the new grid and find the 
values of the field inside the structure of a complex shape. 
The speed of construction of the new matrix depends on the 
dimensions of the structure and the dimensions of the grid. The 
considered method makes it possible to avoid bulky calculation 
of matrix elements. This method is especially efficient for 
calculation of a series of problems for bodies of various 
shapes. Because of the computationally intensive nature of 
the formulated problem, the matrix for the case of a filled 
section has been calculated at the computational cluster of the 
Research Computing Center of the Moscow State University. 





(a) 


Fig 2. The body shape of Q 


Let us present the calculated electric current on a structure 
of a complex shape (Fig. 2).. The results have been obtained 
with the use of a8 x 8 x 8 grid. The incident wave propagates 
along the Ozaxis. The right hand side of the matrix equation 
is E? = Asin (==) Ee N13 | 

The shape of the body is displayed in Fig. 2a, and its 
material the parameter are shown in Fig. 2b. The absolute 
values of the second component of the electric field in the 
first, third, sixth and eighth layers are shown in Figs. 3a— 
3d, respectively. Thus, a subhierarchical method has been 
proposed for solving the problem of diffraction by a dielectric 
body that has an arbitrary shape and is located in a waveguide. 
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(c) (dy 


Fig 3. The modulus of solution of integral equation on the 
body Q 


The method has been implemented on a supercomputer 
calculation, and the are memory required to store the matrix 
is 38 MB for 8x8x8 grid. Program execution time is three 
minutes using 512 processors. We have compared these results 
with an analytic solution to the diffraction problem for (full) 
homogeneous diaphragm in the waveguide. The relative error 
estimation is approximately 1 — 2% for the magnitude of the 
electric field. 


V. CONCLUSION 


The advantages of the method presented in the paper are 
as follows: 1) The integral equation is solved only on the 
bounded body in spite of the fact that the original boundary 
value problem is formulated in the infinite domain; 11) Orders 
of the matrices in the problem are much less than orders of 
similar matrices in the finite difference m ethod or the finite 
elements method; 111) Since coefficients of the matrices can 
be calculated independently, it provides us an opportunity to 
use parallel algorithms to determine the unknown coefficients; 
hence these coefficients can efficiently computed by using 
supercomputers, as has been done in this work. 


This method can be generalized to handle similar 
problems involving waveguides with circular cross-sections. 
We further believe that this approach can be implemented in 
(commercial) software packages to solve a variety of 
diffraction problems, involving inhomoge-neous objects of 
the arbitrary shape. 
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